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Abstract 


A three-stage Runge-Kutta (RK) scheme with multigrid and an implicit preconditioner has been shown to be an 
effective solver for the fluid dynamic equations. This scheme has been applied to both the compressible and 
essentially incompressible Reynolds-averaged Navier-Stokes (RANS) equations using the algebraic turbulence 
model of Baldwin and Lomax (BL). In this paper we focus on the convergence of the RK/implicit scheme when 
the effects of turbulence are represented by either the Spalart-Allmaras model or the Wilcox k-u> model, which 
are frequently used models in practical fluid dynamic applications. Convergence behavior of the scheme with 
these turbulence models and the BL model are directly compared. For this initial investigation we solve the flow 
equations and the partial differential equations of the turbulence models indirectly coupled. With this approach 
we examine the convergence behavior of each system. Both point and line symmetric Gauss-Seidel are considered 
for approximating the inverse of the implicit operator of the flow solver. To solve the turbulence equations we 
use a diagonally dominant alternating direction implicit (DDADI) scheme. Computational results are presented 
for three airfoil flow cases and comparisons are made with experimental data. We demonstrate that the two- 
dimensional RANS equations and transport-type equations for turbulence modeling can be efficiently solved with 
an indirectly coupled algorithm that uses the RK/implicit scheme for the flow equations. 


Contents 


1 Introduction 2 

2 Turbulence Models 3 

2.1 Spalart-Allmaras Model 3 

2.2 Wilcox k-uj Model 4 

3 Numerical Schemes 6 

3.1 RK/implicit Scheme 6 

3.2 DDADI Scheme 8 

4 Computational Results 10 

4.1 Point Symmetric Gauss-Seidel 11 

4.2 Line Symmetric Gauss-Seidel 13 

4.3 Effect of Strain Rate 15 

4.4 Effect of Full Multigrid 16 

5 Concluding Remarks 16 

Appendix 17 


1 



1 Introduction 


Reliable and sufficient convergence for steady-state computations of turbulent flows continues to be a challenge 
in computational fluid dynamics. Here sufficient convergence means that the residuals of the fluid dynamic 
equations and the equation set of a turbulence model are reduced to a level below the truncation error of the 
numerical scheme. In many applications a turbulence model has one or more partial differential equations (PDEs) 
which have a transport form and represent the effects of turbulence on the flow. When solving the transport- 
type equations of turbulence models, either directly or indirectly coupled to the flow equations, the residuals 
are frequently reduced only two orders of magnitude. In addition, the poor convergence of these transport-type 
equations adversely effects the convergence of the flow equations. Of course, when adequate convergence is not 
achieved, there is no assurance that the results obtained represent an acceptable approximation of the solution 
even from an engineering perspective. Thus, there is a strong need for improved numerical methods for not only 
obtaining steady-state solutions but also unsteady solutions when using a dual time-stepping scheme. 

When attempting to develop an improved numerical method for solving the Reynolds-averaged Navier-Stokes 
(RANS) equations a necessary consideration is the coupling of the RANS equations and the equation or equations 
of the turbulence model being applied. If both the fluid dynamic and turbulence equations are directly coupled, 
then the characterization of the discrete system can change. That is, with appropriate discretization the fluid dy- 
namic equations are positive definite (sometimes called a vector positive system [1]), making them amenable to 
relaxation, but the directly coupled system may not be due to the equation set for the turbulence model [2]. The 
numerical stiffness of the entire system is also much higher due to the source terms of the turbulence model. An 
alternative is to use indirect coupling of the two equation sets. Generally, in an iterative solution process with 
this approach the flow variables are updated while the turbulence variables are frozen; and then, the turbulence 
variables are updated while the flow variables are treated as fixed quantities. By indirectly coupling the equations 
we can focus on the specific properties of each equation set to obtain the best possible convergence of the two 
systems of equations. Furthermore, we can identify the essential properties of an algorithm for efficiently solving 
the directly coupled system. There are common design criteria for the algorithms of both equation sets. These 
requirements are as follows: (1) high Courant-Friedrichs-Lewy (CFF) limit, (2) convergence with weak depen- 
dency on mesh density (e.g., computational effort approaching at least 0(N log N)), (3) suitable for stiff discrete 
systems. In addition, for the equation set of the turbulence model there must also be appropriate treatment of any 
source terms so that convergence is not adversely affected. 

A candidate for the flow solver of the loosely coupled system is the RK/implicit scheme. Previously, Rossow [3] 
and Swanson et al. [4] demonstrated that fast convergence can be obtained for both the two-dimensional (2-D) and 
three dimensional (3-D) RANS equations with the RK/implicit scheme and multigrid when using the Baldwin- 
Fomax (BF) algebraic eddy viscosity model [5]. The underlying three-stage RK scheme of this algorithm is 
important for clustering of the eigenvalues associated with the error components of the iterative process. This 
makes the scheme amenable to Krylov subspace iterative methods, which can provide additional convergence 
acceleration. Preconditioning with a fully implicit operator allows a CFF number of 1000. The implicit operator 
can be approximately inverted with either point or line symmetric Gauss-Seidel (SGS). 

The main purpose of this work is to initiate an effort to satisfy the need to significantly augment the effec- 
tiveness (as measured by reliability and efficiency) of algorithms for solving the RANS equations and the PDEs 
of turbulence models. In this preliminary effort we assess the performance of an efficient RANS solver (i.e., 
RK/implicit scheme) when the turbulent viscosity field is generated by solving transport-type equations. Three 
turbulence models are considered. One is the BF model, which is used as a reference model. The other two 
models are the Spalart-Allmaras (SA) and Wilcox k-u> models, which are transport-type equation models that are 
frequently used in solving a variety of fluid dynamics problems. In the first section of this paper we describe these 
two turbulence models indicating the specifics of their implementation. Then the numerical schemes for solving 
the flow and turbulence equations are presented and briefly discussed. For this initial work we use a diagonally 
dominant alternating direction implicit (DDADI) scheme to solve the equation sets of the turbulence models. A 
subiterative procedure is also used to reduce factorization error of the DDADI scheme. In the results section we 
investigate the convergence behavior of the RK/implicit scheme with multigrid when solving the RANS equations 
with the three different turbulence models. The convergence behavior of the scheme with the SA and k-uo models 
is compared with that of the BF model, which was used in the initial development of the numerical scheme. 
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2 Turbulence Models 


The current implementation of the BL model is given in the paper by Swanson and Turkel [6], In the first part 
of this section we present the Spalart-Allmaras model and discuss the particular form of the model used in this 
investigation. Then we describe the Wilcox k-uj model and delineate the specific assumptions associated with 
the form of the model that we have applied. A stress limiter (with parameter) for the specific dissipation rate is 
defined. In a subsequent section on results we discuss the effect of the stress-limiter parameter on the extent of 
flow separation. 


2.1 Spalart-Allmaras Model 

Here we provide a sufficient description of the SA model to allow implementation. A detail discussion explaining 
the modeling of the physical terms in the single transport-type equation is given in the paper by Spalart and 
Allmaras [7]. Let v t be the eddy viscosity, which is defined by 

X 3 

V t = i>fv 1, fvl = 3 r ,3 , x=~, (1) 

X + Wi v 

where v is the kinematic viscosity. The transport-type equation for v given in [7] is written as 


dv 
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dv 
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where t is time, Xj and Uj are Cartesian coordinates and velocity components, respectively, and 


S = S+^f v2 , /»2 = 1 -ttV. (3) 

K z d z 1 + XJv 1 

with S being the magnitude of the vorticity (|fi|), d delineating the distance to the closest wall boundary, and k 
denoting the von Karman constant. The first, second, and third terms on the right-hand side of Eq. 2 represent the 
production, diffusion, and destruction terms, respectively. The last term is a source term, which is defined by 


S=f t iAU 2 , (4) 

where AU is the norm of the difference between the velocity at the transition location and that at a field point 
being considered. The function 


where 


fw — 9 
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g = r + C w2 (r 6 -r), ( 6 ) 

For large values of r the function /„, goes to a constant. Thus, it can be truncated to a value of 10. The function 
f t 2 is defined as 

ft 2 = C t3 exp (— C mX 2 ) O) 

Spalart includes the transition function given by 


ft 1 
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where dt is the distance from the field point to the boundary-layer trip (where trip refers to a known location for 
transition), u>t is the wall vorticity at the trip, 

9t = min [0.1, AU/(u> t Ax t )\. (9) 

and Ax t is the grid spacing along the wall at the trip. 

In the present implementation of the model we do not include the trip function, which is usually neglected 
when applying the model (e.g.. Ref. [8]). In addition, for convenience, after some algebra and rearranging of 
terms, we rewrite Eq. 2 as 


dv dv 
dt Uj dxj 


c b t(i - f t2 )\n\v + {c bl [(i - f t2 )f v 2 + f t2 } K ~ 2 
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( 10 ) 


The constants of the model are as follows: 


C bl = 0.1355, a=\, C b2 = 0.622, k = 0.41, C wl = — + 1 + ° b2 , (11) 

3 K <7 

C w2 = 0.3, C w 3 = 2, C v 1 = 7.1, C*1 = 1, Ct 2 = 2, C t3 = 1.2, C u = 0.5. 

On a solid boundary v = 0. Originally, the free-stream v was set to 1.342 v^, where z/^, is the free-stream 
kinematic viscosity. In order to avoid the possibility of a delayed transition, the free-stream value of v is set to 
3 t'oo, as suggested by Rumsey [9]. 


2.2 Wilcox k-oo Model 


A detailed discussion of Wilcox k-u> model is given in the 3rd edition of the book on turbulence modeling by 
Wilcox [10]. Before writing the two transport equations of this model, we define certain quantities. In general, 
the mean-molecular-stress tensor t, :) and the Reynolds-stress tensor t ? j are given by 


t jj — 2/xSl 


iji 


pTjj 2 fit S^j ^pkSij, 


q Cf _ 1 

ij ~ ij _ 3 dxk ir 


( 12 ) 


where p is density, p is the molecular viscosity, p, t is the turbulent viscosity, k is the turbulence kinetic energy, 
and 8jj is the kronecker delta. The overbar and tilde indicate Reynolds-averaged and mass-averaged quantity, 
respectively. The mean-strain-rate tensor S l3 in Eq. 12 and the mean-rotation tensor f \j are defined as follows: 
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The turbulent viscosity is given by 
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where u is the specific dissipation rate. The equation for the turbulence kinetic energy can be written as 
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The equation for the specific dissipation rate is given by 
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(16) 
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The first, second, and last terms of the right-hand sides of Eqs. 15 and 16 are due to the production, destruction, 
and diffusion of turbulence, respectively. The third term in Eq. 16 is a cross-diffusion term. For the basic imple- 
mentation of the k-u> model we use the the 1988 version [11], following the approach of Ref. [8], This version 
uses the incompressible assumption, which means that the volumetric dilatation term in the total strain-rate tensor 
of Eq. 12 vanishes, and it does not include the cross-diffusion term. In addition, it neglects the term containing k 
in the Reynolds-stress tensor, which is important for high Mach number flows, and it replaces the product of the 
mean-strain-rate tensor and velocity gradient in the turbulence production term by the square of the magnitude of 
the vorticity. 


|ft| 2 = 2 (17) 

where is the zero trace of the mean rotation tensor. When we write it means the scalar (or double dot) 

product of two tensors. Then Eqs. 15 and 16 are replaced by 
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The basic version of the k-co model is modified by using key model changes according to Wilcox [12], For 
additional details see Wilcox [10]. One modification is the introduction of a stress limiter in the definition of co. 
With this change the original co appearing in Eqs. 14, 18, and 19 is replaced in the present form of the k-co model 
by 


lo = max 




C Um = 0.95. 


( 20 ) 


This change is quite important for separated flows. In particular, it reduces the magnitude of the eddy viscosity 
when the production of turbulence energy exceeds the dissipation rate, which allows larger separation bubbles (as 
demonstrated by Huang [13]). The remaining modifications in the model involve the closure coefficients. 

The closure coefficients for the k-co model are as follows: 


/3 = /?o//3, P* = a= \ a * = \' ado= \’ /3 0 = 0.0708, 
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where Prt is the turbulent Prandtl number. The boundary conditions for the model are the following. On a solid 
wall boundary k = k w , which is zero, and lo is determined by 


LO = C0 W 


60/ii 

PiP{di) 2 ’ 


( 22 ) 
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where pi and /j-| are the values of density and molecular viscosity at the center of the mesh cell adjacent to the 
wall. In the free-stream 


k = koo = 9 x 10- 9 i4,, w = Woo = 1 x 10 -6 ( — ) . (23) 

In general, an important consideration when applying turbulence models is the effect of initialization on the 
numerical solution process. For example, when initializing a computation with free-stream conditions, numerical 
difficulties (e.g., convergence stall or even divergence of the iterative procedure for steady-state solutions) can 
occur in regions where uj becomes small, which allows disturbances in the strain rate to result in large values of 
the turbulent viscosity. Menter [14] considers a transport-type equation for the eddy viscosity v t to show that the 
production term becomes very large when uj goes to zero and //, is finite. In order to prevent this from occurring, 
especially during the initial stages of a calculation, Menter suggests applying a constraint on the production term. 
Let P k and D k denote the production and destruction contributions, respectively, in Eq. 20. Then, according to 
Menter a suitable limiter on the turbulence production is given by 

P k = min(P k , 20 D k ). (24) 

We have found that this limiter works well in turbulent airfoil flow applications. 


3 Numerical Schemes 

To solve the RANS equations we use the RK/implicit scheme. Complete details of the scheme are presented in the 
papers of Rossow [3] and Swanson et al. [4]. The SA and k-uj turbulence models require the solution of one and 
two transport-type equations, respectively. In the present work we do not directly couple the solution of the fluid 
dynamic equations with the additional equation or equations of the turbulence models. To solve the transport-type 
equations of the SA and k-uj turbulence models we use a DDADI scheme. In the first part of this section we 
present the essential elements of the RK/implicit scheme. Then we describe the DDADI scheme. 


3.1 RK/implicit Scheme 

We apply a finite-volume approach to discretize the fluid dynamic equations and use the approximate Riemann 
solver of Roe [15] to obtain a second-order discretization of the advection terms. The viscous terms are discretized 
with a second-order central difference approximation. To obtain an explicit update to the solution vector for the 
flow equations we use a three-stage RK scheme. The update for the </-th stage of a RK scheme is given by 

W (9) = W (0) -F <5 W ( «), (25) 


where the change in the solution vector W is 

6W {q) = W (9) - W (0) = -a q ^C W* 9 " 1 *, 


(26) 


and C is the complete difference operator for the system of equations. Here a q is the RK coefficient of the g-th 
stage. At is the time step, and V is the volume of the mesh cell being considered. To extend the support of the 
difference scheme we consider implicit residual smoothing. Applying the smoothing technique of Ref. [16] we 
have the following: 


Lit 5W {q) = SW (q \ 

where L -, is an implicit operator. By approximately inverting the operator L, we obtain 

5W (9) = -a q ^V CW^ = ^ S, 

all faces 


(27) 


(28) 
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where V is a preconditioner defined by the approximate inverse C i 1 , F„ is the normal flux density vector at the 

cell face, and S is the area of the cell face. The change ciW '' q ' replaces the explicit update appearing in Eq. (25). 
Thus, each stage in the RK scheme is preconditioned by an implicit operator. 

A first order upwind approximation based on the Roe scheme is used for the convection derivatives in the 
implicit operator. To derive this operator one treats the spatial discretization terms in the flow equations implicitly 
and applies linearization. For a detailed derivation see Rossow [3], Substituting for the implicit operator in 
Eq. (27), we obtain for the g-th stage of the RK scheme 

W (q) = ~a q ^- Y F^- 1 )5 = R ( «- 1) , (29) 

all faces 

where the matrix A n is the flux Jacobian associated with F„ at a cell face, R ,: 7-1 1 represents the residual function 
for the (g — l)-th stage, and e is an implicit parameter. The parameter e is taken to be 0.7. Analysis and discussion 
concerning choosing e is given in Swanson et al. [4], The matrix A„ can be decomposed into A+ and A”, which 
are associated with the positive and negative eigenvalues of A n and defined by 

A n = 2 (-A-n + |A„|), A n = - (A„ — |A„|). (30) 

If we substitute for A n in Eq. (29) using the definitions of Eq. (30), then the implicit scheme can be written as 


/ + £— Y A n S 

all faces 


I + £ ^~ X A « 5 ' 


all faces 


m!’> = R<y>-e£ -£ A-5W“ s, 


all faces 


(31) 


where the indices (i, j) indicate the cell of interest, and NB refers to all the direct neighbors of the cell being 
considered. Consider traversing the boundary of a mesh cell in a counterclockwise manner, with the positive 
surface normal vector always pointing outward from the cell. Then, as discussed by Rossow [17], the quantity 
A ;; <i W represents the flux density change associated with waves having a negative wave speed (i.e., waves that 
enter the cell (i. j) from outside). Only the neighbor cells NB can contribute to these changes in flux density. 
Similarly, the quantity A+5W represents flux density changes associated with positive wave speeds (i.e., waves 
that leave the cell ( i,j )). These flux density changes are determined only by information from within the cell 

(m)- 

To solve Eq. (31) for the changes in conservative variables hW - ^ , the 4x4 matrix on the left-hand side 
of Eq. (31) must be inverted. It is sufficient to approximate the inverse of the implicit operator. An adequate 
approximate inverse is obtained with two symmetric Gauss-Seidel sweeps. To initialize the iterative process the 
unknowns are set to zero. 

We consider two approaches for approximating the inverse of the implicit operator. One method is to apply 
two sweeps of point symmetric Gauss-Seidel (SGS) relaxation. The other is to apply radial line (or j-line) SGS 
sweeps. We should point out that when referring to this type of SGS in the text the modifiers line and j-line are 
used interchangeably. When computing boundary layer and free shear-layer turbulent flows the computational 
grid is usually strongly clustered to adequately resolve the steep flow-field gradients, which produces a strong 
grid anisotropy and a high degree of stiffness in the discrete flow equations. Line solvers are frequently employed 
with structured grid solvers in the direction of strongest coupling to relieve the stiffness associated with the 
grid anisotropy. Even on unstructured grids line solvers can be applied. By constructing artificial lines in an 
unstructured grid that extend only through the regions of strongest coupling, Mavriplis [18] has demonstrated 
effective reduction of the numerical stiffness. 

With line SGS we apply the implicit operator across the wake line (the line emanating from the trailing edge 
of an airfoil), which we call a fully implicit treatment. It is also possible to lag the boundary condition at the 
wake line rather than extending the line solve across the wake line. This approach can also work well provided 
local relaxation or subiterations are performed to compensate for the lagging of the wake boundary condition. 
Although this may be a more convenient option, from the implementation point of view, for multiblock structured 
grid methods, we use the fully implicit method in this paper to obtain a measure of the best possible convergence 
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with the line SGS relaxation as the smoother for the multigrid method. Furthermore, a variation on this wake 
treatment can be applied with unstructured grid methods by appropriate ordering of the unknowns. 

Due to the upwind approximation used for the implicit operator, the coefficients for the RK scheme are also 
based on an upwind scheme. For computational efficiency we apply a three-stage RK scheme with the coefficients 

[cti, a 2 , < 33 ] = [0.15, 0.4, 1.0] 

from Ref. [19]. In the application of the three-stage RK/implicit scheme as the smoother of a full approximation 
storage (FAS) multigrid method, the CFL number is increased to 1000 after 10 multigrid cycles. During the first 
10 cycles the CFL is 16. A W-type cycle is used to execute the multigrid. Details of the multigrid method are 
given in the paper by Swanson et al. [4]. 

3.2 DDADI Scheme 

For simplicity we describe the DDADI scheme from the perspective of the SA equation. After discretizing Eq. 10, 
we consider the implicit form 


(/ + L x + Ly + S) AF = -R (32) 

where L x and L y are the linear discrete operators for the terms of the transport-type equation, S is the source 
term containing the production and destruction of turbulence contributions, and R is the residual function. The 
operators for the two coordinate directions are as follows: 

L x =Y e l 6 *- (kft + foS x )8 x ] - L v = Y 0 [Sy - (SyP 1 + PiSy)S v ] , (33) 

where S u is a first-order upwind operator for the convection term, 6 is a standard central difference operator, and 
the coefficients [R , fa are defined by the diffusion term of the turbulence model. The parameter 9 indicates the 
temporal accuracy. If 9 = 1/2, then the time derivative is approximated by a central difference, which is second- 
order accurate (i.e., Crank-Nicolson scheme). When (9 = 1 the approximation is a first-order backward difference, 
and we have a fully (an Euler) implicit scheme. The parameter 9 may also be viewed as a measure of implicitness 
with 0 < 9 < 1 and 9 > 1 indicating under-relaxation and over-relaxation, respectively. The source term and the 
residual function are defined by 


5 =—05, R=y K (34 > 

For the advection and diffusion terms of the residual function we use first-order upwind difference and central 
difference approximations, respectively. A first-order approximation of advection terms is frequently applied in 
the implementation of turbulence models in order to promote positivity of the turbulence variables. In general, 
this is not sufficient to ensure positivity, so usually there is also limiting (clipping) of the turbulence quantities 
and/or certain terms (e.g., production term) in the set of turbulence field equations. 

An appropriate linearization of the source term is extremely important to allow the use of large CFL numbers. 
One approach for solving Eq. 32 is to factor the implicit (left-hand side) operator and apply the DDADI scheme. 
Define the diagonal contribution in Eq. 32 as 


D — I + D x + Dy + 5, (35) 

where D x and D y are the diagonal parts of L x and L y , respectively. Then, after factoring out D, we factor the 
resulting operator, obtaining 

D[I + D-\L y - D y )} [I + D~\L X - D x )\ AF = -R. (36) 

The next step is to factor D~ 1 from each term enclosed by brackets, which gives 


[I + Ly + D x + S ] D 1 [I + L x + Dy + S ] AD = -R 


(37) 



To invert this implicit operator, we solve iteratively the one-dimensional systems 


[i + L y + D x + 5] qi — —R, (38) 

[i + L x + D y + S ] <72 = D qi , 

where <72 = AT and 

~n+l =i> n + A ~_ ( 39 ) 

In order to prevent deterioration in the allowable CFL number and damping behavior of the DDADI scheme 
due to the factorization error and possible boundary condition lagging error, we use subiteration. Following 
Klopfer et al. [20] the subiterative process can be written as 


[I + Ly + D x + s ] A<7! = D(q™ - q?) , (40) 

[I + L x + D y + S] A<72 = D(q™ +1 - q™). 

where 

Aq k = C +1 ~€,k = 1, 2; m > 1, ql ? = Av m (41) 

with to denoting the subiteration index. Convergence with iteration and subiteration, Eqs. 38 and 40, can be 
further accelerated by choosing an appropriate implicit parameter 9, just as we do for the fluid dynamic equations. 

The unknown v appears in the denominator of some of the contributions to the source term of the SA model. 
Since the unknown vanishes at the wall, it can happen during the iterative process that at the first point off the wall 
the v becomes zero. This possibility can easily be prevented in the usual way (i.e., take maximum of unknown 
and a very small number). 

For the SA model, four iterations and one subiteration of the DDADI scheme are performed on each stage of 
the RK scheme when on the fine grid of each multigrid cycle for the flow equations. After 15 multigrid cycles the 
number of iterations is reduced to two. In evaluating the fine-grid residuals of the flow equations for restriction, 
we also update the turbulent viscosity. The additional work required for a subiteration is significantly less than 
for an iteration because the implicit coefficients do not have to be recomputed. There is only about a 1% increase 
in computational time per iteration for a subiteration 

As indicated previously, the DDADI scheme just described is also used to solve the two equations of the k-uj 
model. When solving these equations we experience difficulties in simply using a CFL number of 1000 from the 
beginning of the iterative process, as we do when solving the SA model equation. Instead, we increase the CFL 
number from 20 to 1000 over the first 3 multigrid cycles in advancing the solution of the flow equations. For each 
stage of the RK scheme five iterations and one subiteration are performed with the DDADI scheme. 

A desirable property for an effective implicit solver of PDEs is unconditional stability. This property is es- 
pecially important for stiff equations. A linear analysis of an algorithm for solving a discrete system can often 
provide valuable insight into its stability and damping characteristics. MacCormack and Pulliam [21] investi- 
gated the characteristics of the modified approximate factorization (MALk) scheme, with ”k” subiterations per 
time step, by applying Fourier analysis to the Euler equations. The MAF1 scheme is equivalent to the DDADI 
scheme introduced by Bardina and Lombard [22], Their analysis shows that when 9 = 1 the MAF1 scheme is 
conditionally stable. For high CFL numbers instability occurs for small wave number (long wavelength) modes. 
By increasing 9 to 2, MacCormack and Pulliam also show that the DDADI scheme with a single subiteration is 
unconditionally stable. We have observed some different characteristics when solving the PDEs of the SA and 
k-u) turbulence models. In particular, based upon numerical experiments in solving these equations for airfoil 
applications we have demonstrated that the DDADI scheme with one subiteration and 9 = 1 is stable for high 
CFL numbers. With 9 = 2, we need two additional subiterations to maintain stability. 

Figure 1 shows the effect of 9 on the convergence of the DDADI scheme when applying the SA model in a 
transonic airfoil flow computation. Multigrid cycles are used on the horizontal axis for convenience in comparing 


9 




Figure 1: Effect of 9 on the convergence of the DDADI scheme when solving the SA equation (grid: 640 x 128). 


convergence of the flow equations with that of the turbulence equation. For these results the solution to the SA 
equation was updated on each of the three stages of the RK scheme used to solve the flow equations. With each 
update 4 iterations and 3 subiterations per iteration of the DDADI scheme were performed. As revealed in the 
figure, the asymptotic convergence rates are about the same for the three values of 6 being considered. Although 
the residual is reduced slightly more when 9 = 1 and 9 = 2 , the convergence history for 9 = 0.7 exhibits a 
monotonic behavior after the initial iterations. Thus, we generally choose 9 = 0.7. 

When solving the PDEs of the turbulence models along the lines normal to the wake cut line, there are two 
distinct possibilities. One is to treat the wake implicitly, which means to solve for all the unknowns of a model 
along a given radial line in the wake; and thus, no boundary condition is imposed at the wake line. Venkatakrish- 
nan [23] stresses the importance of the fully implicit approach. As indicated previously, an alternative approach 
is to lag the boundary condition at the wake cut line. In general, this method will slow down convergence; how- 
ever, subiteration can be used to eliminate the boundary condition lagging error and convergence slowdown. Here 
again, in the interest of obtaining an estimate of the best convergence, we apply the implicit wake treatment. 

The actual convergence of the DDADI scheme is not independent of the convergence of the flow solver. During 
the course of this investigation we have observed the following convergence behavior. Walsh and Pulliam [24] 
have also reported similar observations. The rate of development of the turbulence field can signficantly affect 
the convergence of the flow solver. Conversely, how well the flow solver converges can have an impact on the 
effectiveness of the scheme for solving the equation set of the turbulence model. Moreover, when the RANS and 
turbulence equations are being solved in a loosely coupled manner, an essential requirement for an effective total 
algorithm is that the numerical solution vector of each equation set exhibits a similar evolution rate. Since the 
present algorithm includes an effective flow solver, we have provided an amenable setting for the DDADI scheme 
and the overall algorithm. 

4 Computational Results 

A series of computations for turbulent, viscous flow over the RAE 2822 airfoil were performed to evaluate the 
convergence behavior of the RK/implicit scheme when applying the BL, SA, and k-to turbulence models. In solv- 
ing the RANS equations we used structured meshes with a C-type topology. We considered two mesh densities, 
one with 320 cells around the airfoil (256 cells on the airfoil) and 64 cells in the normal direction, and the other 
with twice as many cells in each coordinate direction. For the 640 x 128 grid the normal mesh spacing at the 
surface is approximately 6 x 10 , and the maximum surface cell aspect ratio is about 2000. The airfoil solu- 
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Table 1 : Flow conditions for RAE 2822 airfoil. 


Cases 

Mx, 

a (deg.) 

R&c 

%tr / C 

Case 1 

0.676 

1.93 

5.7 x 10 6 

0.11 

Case 9 

0.730 

2.79 

6.5 x 10 6 

0.03 

Case 10 

0.750 

2.81 

6.2 x 10 6 

0.03 


tions were calculated with the Case 1, Case 9 and Case 10 flow conditions (see Table 1) from the experimental 
investigation of Cook, McDonald and Firmin [25]. For Case 1 the flow is primarily subsonic with a relatively 
small region of supersonic flow. The other two cases are transonic. In Case 9 there is a fairly strong shock wave 
occurring on the upper surface of the airfoil, where as in Case 10 there is a sufficiently strong shock on the upper 
surface to cause significant flow separation behind the shock. 

In all the applications the same boundary conditions were imposed for the fluid dynamic equations. On the 
surface the no-slip condition was applied. At the outer boundary Riemann invariants were used. A far-field vortex 
effect was included to specify the velocity for an inflow condition at the outer boundary. A detailed discussion 
of the boundary conditions is given in Ref. [6], The calculations were started on the solution grid with the initial 
solution given as the free-stream conditions. When a full multigrid process was applied, the initial solution on a 
given level of refinement was obtained from a coarser level. 

When solving the fluid dynamic equations with the RK/implicit scheme both point and line SGS were con- 
sidered to approximate the inverse of the implicit operator. With the multigrid method four grid levels were used. 
Unless indicated otherwise, full multigrid (FMG) was not applied when solving the flow equations, making a 
more demanding challenge for the initial part of the iterative process. All computations were performed on a Dell 
XPS 630 computer, which has two quad-core processors at 2.66 GHz. 


4.1 Point Symmetric Gauss-Seidel 

In this section we investigate the convergence behavior of the RK/implicit scheme with point SGS(2), where the 
2 refers to the number of symmetric sweeps, for approximating the inverse of the implicit operator. We also 
examine the effect of limiting the aspect ratio of the mesh cells in the wake of the airfoil flow. The L 2 norm of the 
residual of the continuity equation is used as a measure of convergence for the flow equations. For these results 
the calculations were terminated after either 125 multigrid cycles or a 13 orders of magnitude reduction in the 
residual. 

We first consider the convergence of the flow equations with the RK/implicit scheme when applying the BL 
model. Figure 2 shows convergence histories of the residual and lift coefficient (Cl) for Case 1 on both the 
320 x 64 and 640 x 128 grids. On the 320 x 64 grid the scheme converges at a rate of 0.714. There is a some 
slowdown in convergence on the 640 x 128, where the rate is 0.742. With the turbulent viscosity field produced 
by the BL model there is clearly no reduction in the asymptotic convergence rate (i.e., rate beyond approximately 
8 orders of magnitude reduction in the residual) on a given grid. As we will see next, the turbulent viscosity fields 
generated with the other turbulence models do cause a slower asymptotic convergence rate when using point 
SGS(2). 

Figure 3 shows the convergence of the solvers of the flow equations and the SA transport-type equation for 
Case 1 on the 320 x 64 grid. The CFL number for both the flow equations and the turbulence equation was 1000. 
As seen in Figure 3(a), there is significant slowdown in convergence after about 40 multigrid cycles (residual 
reduced 8 orders) in the absence of no limiting of the mesh cell aspect ratio (AR) in the wake. By limiting the 
wake aspect ratio to 1000 the residual is reduced two additional orders before a slowdown occurs. In the limiting 
the normal grid spacing in the airfoil wake is allowed to increase with distance downstream of the airfoil trailing 
edge. With an AR limit of 100 there is no deterioration in the asymptotic convergence rate. The convergence of 
the residual of the SA equation, which is displayed in Figure 3(b), responds in a similar way to the AR limiting. 

Figure 4 shows the convergence histories on the 320 x 64 grid when applying the k-ui turbulence model. For 
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Figure 2: Convergence histories of Case 1 using point SGS with RK/implicit scheme and the BL model (grids: 
320 x 64, 640 x 128). 


the flow equations we again observe the slower asymptotic convergence rate when using the point SGS(2) scheme 
in the multigrid smoother. This slowdown is eliminated for both the flow and turbulence equations by limiting the 
wake AR to 1000. These results also seem to indicate that the k-oj model is less sensitive than the SA model to 
AR in the wake. 

Clearly, for the two models the residual has been reduced well below the truncation error of the numerical 
discretization before the slower residual reduction rate occurs. This means that we have achieved converged 
solutions well before the onset of slowdown, and that further reduction in the residual does not improve the 
solutions. However, in more complex flow problems the source of this slowdown may prevent convergence. For 
this reason it is important to understand the cause of the slower asymptotic rate of convergence. 

We now address the issue of why there is a slowdown in the asymptotic convergence rate when applying either 
the S A or the k-ui models and using point SGS for inverting the implicit operator in the RK/implicit scheme. Since 
there is no apparent slowdown indicated by the asymptotic convergence rate when the BL model and point SGS are 
employed, the mesh AR in the wake does not appear to be the primary source of the slower asymptotic rate when 
applying the SA and k-u> models. Another possible cause of the reduction in convergence rate is the difference in 
the turbulent viscosity fields in the wake. Thus, we consider a characterization of turbulent viscosity for all the 
turbulence models in the airfoil wake. 

An important factor for the different turbulence models is how well they represent the physics of the airfoil near 
wake flow. The manner in which the models behave in the near wake can also significantly effect the convergence 
behavior of both the fluid dynamic and turbulence equations. Here the issue is not only the magnitude of the 
turbulent viscosity (jit) but also how it varies as the flow proceeds from a wall-bounded shear layer to a free 
shear layer. In Figs. 5(a) - Figs. 5(c) we present the ji t contours in the near wake region from solutions on the 
320 x 64 grid obtained with the three turbulence models. The jx t . is scaled by the free-stream molecular viscosity 
(jj, oo). These results indicate that in the near wake the streamwise extent of higher levels of //, is much greater 
for the SA and k-u models than for the BL model. Also, with the k-u model, we observe unexpected rapid 
initial spreading of the ji t field in the initial wake, resulting in a much greater lateral extent. Figure 5(d) shows 
that mesh refinement (640 x 128 grid) produces a significant reduction of the ji t field in the transverse direction, 
suggesting that resolution plays a key role in the development of the near wake flow. Some further discussion of 
the initial wake behavior is given in the Appendix. Turbulent viscosity profiles at four streamwise wake locations 
are displayed for each of the turbulence models in Fig. 6. For the BL model the peak value of \i t occurs at the 
trailing edge, and a minimum value occurs at approximately the two chord location. From the complete turbulent 
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Table 2: Computational performance for Case 1 using different turbulence models. LSGS(2) used in solving flow 
equations and DDADI used to solve turbulence equations. 


Model 

Grid 

CPU Time (sec.) 

MG Cycles 

Convergence Rate 

BL 

320 x 64 

79 

79 

0.683 

SA 

320 x 64 

100 

79 

0.684 

k-u 

320 x 64 

125 

80 

0.686 

BL 

640 x 128 

365 

87 

0.709 

SA 

640 x 128 

544 

87 

0.708 

k-u 

640 x 128 

673 

88 

0.710 


viscosity field we observe that a maximum fj, t value of about 400 persists into the far field. In contrast, both the 
SA and k-u models have their peak value of // t occurring at around the two chord location, and then the peak fj, t 
decreases as the flow proceeds downstream. The maximum ji t of these two models is roughly twice as large as 
that of the BL model. Figure 6 also shows that the fi t values in the vicinity of the y/c = 0 line (i.e., the wake cut) 
are much higher for the SA and k-u models than for the BL model. These higher levels of /q result in a much 
stronger coupling of the discrete equations through the diffusion terms. In addition, there is an augmentation in 
the influence of the source terms in the turbulence models. Both of these factors contribute to the increased need 
for implicit treatment across the cut line of the wake. 


4.2 Line Symmetric Gauss-Seidel 

In the next set of results we show the benefit of using j-line symmetric Gauss-Seidel (LSGS) to approximate the 
inverse of the implicit operator in the RK/implicit smoother. Figure 7 compares the convergence histories of the 
flow equations for Case 1 on the two grids when applying the three different turbulence models. With each model 
nearly the same convergence rate is obtained on 320 x 64 grid. Moreover, almost the same rate is attained for each 
model on the 640 x 128 grid. The convergence rates are given in Table 2. Figure 8 displays the convergence plots 
for the SA and k-ui models. For both of these models there is some slowdown in the convergence rate due to mesh 
refinement. In general, slower rates are expected since there is no component in the algorithm for the turbulence 
equations, such as multigrid, to reduce this effect. We need to keep in mind that if multigrid is introduced the 
emphasis is placed on the role of the scheme as a high-frequency smoother rather than a solver. As indicated in 
the papers by MacCormack and Pulliam [21] and Pulliam et al. [26], the DDADI scheme with subiteration has 
potential as an effective smoother. 

The surface pressure and skin-friction distributions computed with the three turbulence models for Case 1 are 
compared in Figs. 9 and 10. The surface pressures are quite similar and exhibit good agreement with the data. 
The skin-friction coefficient C / is the wall shear stress divided by the dynamic pressure based on boundary-layer 
edge values. In Table 3 we give the computed lift and drag coefficients using the three models. 

For Case 9 the convergence plots for the flow and turbulence equations are presented in Figs. 11 and 12. In 
this transonic case we observe somewhat similar convergence histories for the BL and SA turbulence models 
on both meshes. The convergence rate with these models is about 0.7 (see Table 4). With the k-u model the 
convergence behavior is similar to that of the other two models on the 320 x 64 grid, but there is a brief slowdown 
in the convergence at around 30 cycles on the 640 x 128 grid. However, there is a fairly quick recovery in the 
convergence, and the asymptotic rate is similar to that of the other two models. This slowdown appears to be a 
consequence of the stress limiter of the model, which will be confirmed later. The convergence plots for the SA 
model are shown in Fig. 12(a). On the two grids the convergence rates are roughly the same after 40 cycles for 
the flow equations, but then there is an unexpected more rapid residual reduction on the finer grid. In addition, 
we do not see the oscillatory behavior that occurs for Case 1 . For the k-u model the convergence rates on the two 
grids are about the same. Although we do not observe a slower convergence with mesh refinement, one should 
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Table 3: Computed lift and drag coefficients for RAE 2822 airfoil. Numerical dissipation from Roe scheme. Weak 
grid clustering in neighborhood of shock wave. 


Cases 

Model 

Grid 

C L 

c D 

{Cd) p 

{C D )v 

Case 1 

BE 

320 x 64 

0.6096 

0.008280 

0.002490 

0.005790 

Case 1 

BE 

640 x 128 

0.6096 

0.008182 

0.002380 

0.005802 

Case 1 

SA 

320 x 64 

0.5913 

0.008458 

0.002617 

0.005841 

Case 1 

SA 

640 x 128 

0.5902 

0.008272 

0.002486 

0.005785 

Case 1 

k-u> 

320 x 64 

0.6027 

0.008183 

0.002410 

0.005773 

Case 1 

k-uj 

640 x 128 

0.5933 

0.008188 

0.002400 

0.005788 

Case 9 

BE 

320 x 64 

0.8541 

0.01783 

0.01230 

0.005531 

Case 9 

BL 

640 x 128 

0.8613 

0.01804 

0.01246 

0.005585 

Case 9 

SA 

320 x 64 

0.8186 

0.01671 

0.01118 

0.005537 

Case 9 

SA 

640 x 128 

0.8220 

0.01662 

0.01111 

0.005507 

Case 9 

k-LO 

320 x 64 

0.8333 

0.01702 

0.01155 

0.005477 

Case 9 

k-io 

640 x 128 

0.8178 

0.01644 

0.01090 

0.005545 

Case 10 

BL 

320 x 64 

0.8523 

0.02910 

0.02366 

0.005435 

Case 10 

BL 

640 x 128 

0.8781 

0.03067 

0.02512 

0.005553 

Case 10 

SA 

320 x 64 

0.7775 

0.02584 

0.02052 

0.005317 

Case 10 

SA 

640 x 128 

0.7850 

0.02598 

0.02069 

0.005290 

Case 10 

k-LO 

320 x 64 

0.7848 

0.02581 

0.02056 

0.005256 

Case 10 

k-LO 

640 x 128 

0.7672 

0.02484 

0.01951 

0.005331 


anticipate that on finer meshes a convergence slowdown will occur. Again this effect can possibly be limited by 
introducing multigrid into the algorithm for the turbulence equations. Furthermore, with the multigrid it may be 
possible to reduce the number of applications per RK stage of the DDADI scheme or whatever solver is being 
used for the PDEs of the turbulence models. In Figs. 13 and 14 we present the computed surface pressure and 
skin-friction distributions. The primary differences in the results with the three models is in the vicinity of the 
upper surface shock. For the SA and k-ui models the surface pressures are nearly the same. The variation in 
surface skin-friction is also quite similar, with the largest differences occurring in the region downstream of the 
shock. With the BE model the shock location is farther downstream, resulting in higher lift and drag coefficients 
(see Table 3). 

The final set of results is for Case 10. In Figs. 15 and 16 the convergence plots for the flow and turbulence 
equations are presented. This transonic case is more difficult due to the stronger shock, which induces flow 
separation. The flow solver converges faster (rate of 0.72 on the 640 x 128 grid) when using the SA model. 
With the BE model the asymptotic convergence rate is slower. A comparison of the convergence rates is given in 
Table 5. From the convergence plots for the SA model, which are displayed in Fig. 16(a), we see that the residuals 
are reduced at least 9 orders of magnitude, and the rate of convergence on the 640 x 128 grid is somewhat faster. 
For the k-u> model the DDADI scheme converges roughly the same for both the turbulence kinetic energy and the 
specific dissipation rate equations on both grids. In Figs. 17 and 18 we present the computed surface pressure and 
skin-friction distributions. The surface pressures predicted with the BE model show the shock location to be much 
farther downstream than the experiment does. Clearly, there is much better agreement with the data when applying 
the SA model and k-io models. According to Wilcox [12], the k-u model with the stress limiter gives essentially 
the same results for these types of flows as the shear-stress transport (SST) model of Menter [14] does. This can 
be verified for Case 10 by comparing the present result with that calculated by Menter and Rumsey [27] using 
the SST model. As expected the principal difference in the surface skin-friction values computed with the three 
turbulence models is due to the prediction of the shock location. Even though there is a relatively small difference 
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Table 4: Computational performance for Case 9 using different turbulence models. LSGS(2) used in solving flow 
equations and DDADI used to solve turbulence equations. 


Model 

Grid 

CPU Time (sec.) 

MG Cycles 

Convergence Rate 

BL 

320 x 64 

90 

91 

0.720 

SA 

320 x 64 

109 

87 

0.708 

k-LO 

320 x 64 

142 

92 

0.722 

BL 

640 x 128 

343 

81 

0.690 

SA 

640 x 128 

504 

80 

0.688 

k-LO 

640 x 128 

685 

89 

0.714 


Table 5: Computational performance for Case 10 using different turbulence models. LSGS(2) used in solving 
flow equations and DDADI used to solve turbulence equations. 


Model 

Grid 

CPU Time (sec.) 

MG Cycles 

Convergence Rate 

BL 

320 x 64 

105 

105 

0.752 

SA 

320 x 64 

126 

99 

0.739 

k-LO 

320 x 64 

178 

115 

0.770 

BL 

640 x 128 

399 

95 

0.729 

SA 

640 x 128 

587 

93 

0.724 

k-LO 

640 x 128 

888 

116 

0.772 


in the shock position obtained with the SA and k-ui models, the behavior of the flow with the two models is 
different downstream of the shock. The SA model indicates either flow separation or near separation extending 
to the airfoil trailing edge, but the k-to model indicates separated flow only from x/c = 0.55 to x/c = 0.70. 
Moreover, the attached flow solution with the k-ui model appears to be in better agreement with the experiment. 
In Table 3 the computed lift and drag coefficients using the different models are given. 

The stress limiter of the k-u> model has a significant effect on not only the representation of the physics for 
this strong transonic case but also the convergence of both the flow and turbulence equations. Figure 19 shows the 
influence of the stress limiter parameter Cu m on convergence. The Cu m is varied between 0.0 and 0.95, which is 
the value used in the previous results. For this range of Cu rn the number of cycles required to reduce the residual 
of the flow equations 13 orders increases by 35%. The residuals of the turbulence equations are two orders higher 
when Cum = 0.95. It may be possible to improve the convergence through the use of a smoother switching from 
the unlimited specific dissipation rate to the limited one (see Eq. 20). Figure 20 displays the strong effect on the 
shock location with the variation of the parameter Cu m . 


4.3 Effect of Strain Rate 

In all of the previous results we have used the magnitude of the vorticity in defining the turbulence production 
(I\) terms and the stress limiter of the specific dissipation rate in the k-co model (see Eqs. 18, 19, and 20). 
Since Wilcox [10] uses the strain rate S tJ to determine rather than the vorticity magnitude in his version of 
the k-uj model, the question arises as to the effect of this on not only the results but also the convergence of the 
RK/implicit scheme. To address this question we computed the solution of Case 10 with the strain rate. Figures 21 
and 22 show the effect of the two different ways of computing the Pi- terms. In Figs. 21(a) and 21(b) we see that 
there is little difference in the convergence histories for the flow and turbulence equations when using strain rate. 
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Table 6: Comparison of computed lift and drag coefficients for Case 9 on 640 x 128 grid using FMG. For the 
reference result with SA and k-u > models the residual was reduced 13 orders of magnitude. 


Model 

MG Cycles 

c L 

c D 

SA 

63 

0.8220 

0.01662 

SA 

10 

0.8220 

0.01662 

SA 

5 

0.8213 

0.01660 

k-u 

69 

0.8178 

0.01644 

k-u 

10 

0.8179 

0.01644 

k-u 

5 

0.8191 

0.01651 


The comparisons of the surface pressure and skin-friction distributions in Fig. 22 reveal generally only slight 
differences due to the method of calculating the T\ terms. There is a small downstream shift in the indicated 
transition location when using strain rate. Although not shown, the solutions obtained for Cases 1 and 9 also show 
little variation due to the two methods for calculating . We should point out that the Cu m in the stress limiter of 
the Wilcox model is 0.875. With the present implementation of the k-u model, we found that this value resulted 
in the computed shock location being downstream of the measured data. By using Cu m = 0.96, nearly the same 
value used when computing with vorticity magnitude, we obtained good agreement with the data. There are 
other differences between the current implementation of the model and the 2006 Wilcox implementation that may 
possibly explain the different value for the parameter Cu rn (see Section 2.2). 

4.4 Effect of Full Multigrid 

Convergence of the solution, to the approximate level of the truncation error, can be accelerated by implementing 
full multigrid (FMG). We now consider the performance of the RK/implicit scheme with FMG when applying 
the SA and k-u models. The residual and lift coefficient (Cl) convergence histories for the 3-stage RK/implicit 
scheme when applying the SA model and using FMG are shown in Fig. 23. The calculation was done for Case 
9 with the 640 x 128 grid using 3 levels of refinement, which contain 3, 4, and 5 grids. On each level multigrid 
was executed until either 100 cycles were completed or the residual was reduced 13 orders of magnitude. The 
CFL number for the flow equations was increased to 1000 over the first 10 cycles of the first refinement level 
and maintained on all the remaining levels. As seen in Table 6, with only 5 cycles on the final level the Cl 
error relative to the converged value is less than 0.1%. Agreement with the converged Cl and C /> is obtained to 4 
significant figures in 10 cycles. As evident in Fig. 24 the surface pressure and skin-friction distributions at 5 cycles 
is nearly identical to the corresponding one with the residual reduced 13 orders of magnitude. In Figs. 25 and 26 
we present the convergence histories and solution comparisons when applying the k-u model. After 5 cycles both 
the Cl and Cd are within 0.5% of the converged value. Here again to plotting accuracy we see essentially the 
same surface pressure and skin-friction distributions for the 5 cycle and converged solutions. 


5 Concluding Remarks 

In this work we have investigated the performance of the RK/implicit algorithm for solving the RANS equations 
when different types of turbulence models are applied. Airfoil flows have been solved with three turbulence 
models: (1) the algebraic model of Baldwin and Lomax, (2) the one-equation model of Spalart and Allmaras, 
(3) the two-equation (k-u) model of Wilcox. The fluid dynamic and turbulence equations have been solved in 
a loosely coupled manner. This approach has allowed examination of the convergence properties of each set 
of equations. The fluid dynamic equations have been solved with the RK/implicit scheme and multigrid. Both 
point and line solves have been considered to approximate the inverse of the implicit operator. Line solves are in 
general necessary to appropriately treat the stronger coupling of the flow equations in an airfoil wake when the 
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eddy viscosity field is determined with either the SA or k-cu models. The equation set for each turbulence model 
has been solved only on the solution grid with a DDADI scheme. For both the fluid dynamic and turbulence 
equations a CFL of 1000 has been used. 

To evaluate each part of the numerical algorithm for solving the loosely coupled system of equations, we have 
used three of the standard RAE airfoil test cases (Cases 1, 9, and 10). For all cases computational results have 
been obtained on a standard grid (320 x 64 cells) and on a fine grid (640 x 128 cells) in order to indicate the effect 
of mesh refinement on convergence. We have investigated not only the convergence behavior of the algorithm but 
also the effects of the turbulence modeling on the solutions. In all cases the flow equations converged as fast, if 
not faster, when using the SA and k-ui models instead of the BL model. Furthermore, these equations converged 
at similar rates on both the standard and fine grids. For the equation sets of the SA and k-u) models, we have 
shown that the residuals can be reduced about 10 orders of magnitude when the residuals for the flow equations 
are reduced 13 orders. With the DDADI scheme we have not observed a strong effect on the rate of convergence 
due to mesh refinement. 

With all turbulence models reasonable agreement has been obtained with the experimental surface pressure 
and skin-friction distributions for Cases 1 and 9. For Case 10, where there is a sufficiently strong shock to induce 
significant separation, only the SA and k-u> models exhibit good agreement with the measured data. The k-oj 
model requires the introduction of a stress limiter for the specific dissipation rate. With the stress limiter the k-oj 
model can produce essentially the same results obtained with Menter’s SST model. Without this limiter the shock 
location is predicted too far downstream. 

In applying the RK/implicit scheme with j-line SGS a convergence rate of approximately 0.7 has been obtained 
for both Case 1 and Case 9 on the standard and fine grids with all turbulence models. The convergence rates of 
the DDADI scheme on each grid are similar for a given turbulence model. For Case 10 the flow solver has 
converged at rates of about 0.74 and 0.73 with the BL and SA models, respectively. With the k-cu models there is 
some deterioration in the average convergence rate of the scheme (sa 0.77). Yet, the asymptotic rate is about the 
same as obtained with the other models. When solving the PDEs of the SA and k-u models the DDADI scheme 
converges roughly at the same rate as for Case 9. As part of the convergence investigation we have also considered 
the effect of FMG. With only 5 multigrid cycles on the finest grid we have essentially obtained the solution to 
plotting accuracy for Case 9 with both the SA and k-u models. 

At this point we have demonstrated that the 2-D RANS equations and transport-type equations can be solved 
efficiently in a loosely coupled manner. In the continuation of this work we plan to consider several enhancements 
of the current algorithm to improve efficiency and robustness. The DDADI scheme with subiteration was chosen 
in this initial effort to solve the PDE equation set of a turbulence model because it is an effective solver, and it 
allowed a reference efficiency to be established. With any scheme for solving the equation set of a turbulence 
model there will be slowdown in convergence on sufficiently fine meshes. To reduce any convergence effects 
due to mesh density we will include multigrid in the scheme for solving the PDEs of the turbulence models. By 
introducing multigrid the role of what was the solver will be changed to that of a high-frequency smoother. Various 
smoothers will be considered, including the RK/implicit scheme. An acceptable smoother must be effective in 
three dimensions. Ultimately, the objective is to develop both loosely coupled and fully coupled algorithms that 
are reliable and efficient in both two and three dimensions. 

So far we have not directly addressed the issue of robustness when solving the turbulence equations. A 
powerful technique for increasing reliability of the complete algorithm is the Krylov subspace iterative method. 
Frequently, there are only a few error modes that inhibit or stall the convergence of a scheme. Since the Krylov 
method works on the dominant error modes, it can effectively annihilate these troublesome error components. 
The effects of a Krylov method on the performance of the scheme will also be investigated. 

Appendix 

In section 4.1 we showed that the rapid increase in the lateral extent of the near wake /r t field determined with 
the k-u model (see Fig. 5) can be significantly affected by increasing the mesh resolution. Further investigation 
has indicated that accuracy is a controlling factor for obtaining reasonable physical behavior in the near wake 
development of the fi t field. In this investigation we also determined that inclusion of the cross-diffusion term 
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produced only a small effect on the /i t behavior in the initial wake region. Figure 27(a) displays the trailing edge 
region of the 320 x 64 grid used in the previous calculations. For this grid, which was generated by a hyperbolic 
method, there is a spreading out of the radial lines. By choosing different grid generation parameters we generated 
a new grid with a dramatically reduced spreading effect, as shown in Fig. 27(b). We designate the original and 
new grids as the reference and modified 320 x 64 grids, respectively. The influence of the modified grid on the 
fit variation is displayed in Fig. 28. Clearly, there is a significant reduction in the lateral extent of the non-zero / x t 
in the near wake region. Comparing the profiles of Fig. 28 with those of Fig. 6(c) and Fig. 6(d) we see that the Ht 
profiles on a given grid density are actually quite similar, with some differences due to lateral extent. 

Another factor that can effect the near wake behavior is the order of approximation of the convection terms 
in the model. White [28] has computed solutions for Case 1 with both a first-order and a second-order upwind 
approximation of the convection terms in essentially the 2006 Wilcox k-u> model [10], Fie solved the fully cou- 
pled RANS and turbulence model equations with a second-order, finite-volume algorithm [29]. On the present 
reference 640 x 128 grid White obtained very similar //, contours to those of Fig. 5(d) with the first-order approx- 
imation of the convection terms. Using a second-order upwind approximation he obtained similar //, contours 
to those of Fig. 28(c). These results suggests that the convection terms play a much more important role in a 
free-shear flow with the k-w model than the SA model. Consequently, the k-uj model is much more sensitive to 
the accuracy of the convection terms in the wake flow than the SA model. 
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(a) (b) 

Figure 3: Effect of wake aspect ratio on convergence of Case 1 with SA model on 320 x 64 grid, (a) Flow 
equations, (b) SA equation. 




(a) (b) 

Figure 4: Effect of wake aspect ratio on convergence of Case 1 with k-ui model on 320 x 64 grid, (a) Flow 
equations, (b) k-ui equations. 
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Figure 5: Contours of turbulent viscosity near trailing edge for Case 1, computed with BL, SA, and k-u> models. 
On 320 x 64 grid: (a) BL model, (b) SA model, (c) k-u> model. Effect of grid refinement: (d) k-u> model, 640 x 128 
grid. 
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(a) 



(b) 




(c) (d) 


Figure 6: Turbulent viscosity profiles in the airfoil wake for Case 1, computed with BL, SA, and k-ui models. On 
320 x 64 grid: (a) BL model, (b) SA model, (c) k-u model. Effect of grid refinement: fd) k-ui model, 640 x 128 
grid. 
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Figure 7: Convergence histories for Case 1 on two grids with BL, SA, and k-u> models. Flow equations solved 
with line SGS(2) and transport-type equations solved with DDADI. (a) 320 x 64 grid, (b) 640 x 128 grid. 



Figure 8: Convergence histories of SA and k-ui models for Case 1 on two grids. Transport-type equations solved 
with DDADI. (a) SA model, (b) k-u> model. 
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(a) (b) 

Figure 11: Convergence histories for Case 9 on two grids with BL, SA, and k-ui models. Flow equations solved 
with line SGS(2) and transport-type equations solved with DDADI. (a) 320 x 64 grid, (b) 640 x 128 grid. 




(a) (b) 

Figure 12: Convergence histories of SA and k-u> models for Case 9 on two grids. Transport-type equations solved 
with DDADI. (a) SA model, (b) k-u model. 
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Figure 13: Surface pressure distributions for Case 9 on two grids with BL, SA, and k-u> models, (a) 320 x 64 
grid, (b) 640 x 128 grid. 
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Figure 14: Surface skin-friction distributions for Case 9 on two grids with BL, SA, and k-uj models, (a) 320 x 64 
grid, (b) 640 x 128 grid. 









(a) (b) 

Figure 15: Convergence histories for Case 10 on two grids with BL, SA, and k-ui models. Flow equations solved 
with line SGS(2) and transport-type equations solved with DDADI. (a) 320 x 64 grid, (b) 640 x 128 grid. 




(a) 


(b) 


Figure 16: Convergence histories of SA and k-w models for Case 10 on two grids. Transport-type equations 
solved with DDADI. (a) SA model, (b) k-u model. 
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(a) (b) 

Figure 19: Effect of stress limiter of k-u> model on convergence histories of fluid and turbulence equations (Case 
10, 320 x 64 grid). Flow equations solved with line SGS(2) and k-oj model equations solved with DDADI(5,1). 
(a) Flow equations, (b) turbulence equations. 




(a) (b) 


Figure 20: Effect of stress limiter of k-ui model on surface pressure and skin-friction distributions for Case 10 
(320 x 64 grid), (a) Surface pressures, (b) Surface skin-friction. 
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Figure 21: Effect of computing turbulence production (Pk) terms of k-u> model with vorticity or strain rate. 
Convergence histories of fluid and turbulence equations (Case 10, 640 x 128 grid). Flow equations solved with 
j-line SGS(2) and k-u model equations solved with DDADI(5,1). (a) Flow equations, (b) turbulence equations. 
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Figure 22: Effect of two ways of computing Pk terms of k-ui model on surface pressure and skin-friction distri- 
butions for Case 10 (320 x 64 grid), (a) Surface pressures, (b) Surface skin-friction. 
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Figure 23: Convergence history of flow equations with FMG (3 levels of refinement) for Case 9. Flow equations 
solved with j-line SGS(2) and SA model equation solved with DDADI. 
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Figure 24: Comparison of surface pressure and skin-friction distributions at 5 cycles and 63 cycles (residual 
reduced 13 orders) computed for Case 9 with SA model. FMG used to solve flow equations, (a) Surface pressures, 
(b) Surface skin-friction. 
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Figure 25: Convergence history of flow equations with FMG (3 levels of refinement) for Case 9. Flow equations 
solved with j-line SGS(2) and k-w model equations solved with DDADI. 
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Figure 26: Comparison of surface pressure and skin-friction distributions at 5 cycles and 69 cycles (residual 
reduced 13 orders) computed for Case 9 with k-to model. FMG used to solve flow equations, (a) Surface pressures, 
(b) Surface skin-friction. 
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Figure 27: Airfoil trailing edge region of two different grids with 320 x 64 cells. 
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(b) Modified grid. 
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Case 1 , Turbulent Viscosity Contours, k-co Model 
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Case 1 , Turbulent Viscosity Contours, k-co Model 
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Figure 28: Effect of modified grid on turbulent viscosity contours and profiles near the trailing edge for Case 1, 
computed with the k-co model. On 320 x 64 grid: (a) contours, (b) profiles. On 640 x 128 grid: (c) contours, fd) 
profiles. 
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